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Abstract. We seek solutions u £ R n to the semilinear elliptic partial difference equation —Lu + 
f s (u) = 0, where L is the matrix corresponding to the Laplacian operator on a graph G and 
f s is a one-parameter family of nonlinear functions. This article combines the ideas introduced 
by the authors in two papers: a) Nonlinear Elliptic Partial Difference Equations on Graphs (J. 
Experimental Mathematics, 2006), which introduces analytical and numerical techniques for solving 
such equations, and b) Symmetry and Automated Branch Following for a Semilinear Elliptic PDE 
on a Fractal Region (SIAM J. of Dynamical Systems, 2006), wherein we present some of our recent 
advances concerning symmetry, bifurcation, and automation for PDE. 

We apply the symmetry analysis found in the SIAM paper to arbitrary graphs in order to obtain 
better initial guesses for Newton's method, create informative graphics, and better understand the 
role of symmetry in the underlying variational structure. We use two modified implementations 
of the gradient Newton- Galerkin algorithm (GNGA, Neuberger and Swift) to follow bifurcation 
branches in a robust way. By handling difficulties that arise when encountering accidental degen- 
eracies and higher-dimensional critical eigenspaces, we can find many solutions of many symme- 
try types to the discrete nonlinear system. We present a selection of experimental results which 
demonstrate our algorithm's capability to automatically generate bifurcation diagrams and solu- 
tion graphics starting with only an edgelist of a graph. We highlight interesting symmetry and 
variational phenomena. 



This paper considers nonlinear partial difference equations (PdE) on graphs. In particular, we 
automate the bifurcation analysis and branch following required for finding solutions u £ R n to the 
discrete nonlinear system 



Here, L is the matrix corresponding to the Laplacian operator on a simple connected graph G and 
f s : R — > R satisfies / s (0) = and f' s (0) = s. The nonlinear term f s (u) £ R n is defined as a 
composition, that is, (f s (u))i = f s (ui). The real number s is treated as a bifurcation parameter. 
The existence of the trivial solution u = £ R n is clear for all s G R, since / s (0) = 0. By finding 
and following new, bifurcating branches of (generally) lesser symmetry we are able to find, within 
reason, any solution that is connected by branches to the trivial branch. Our code works for a 
wide range of nonlinearities. It does not require that f s is odd, that the nonlinearity is superlinear 
[U [6], nor that f s has the form f s (t) = st + H(t). That being said, in this paper we choose f s to 
be the family of odd and superlinear functions defined by f s (t) = st + i 3 except when otherwise 
specified. Our ultimate goal is to automate the process of accurately approximating all solutions to 
Equation (pQ) given only the edgelist for a given graph G, and then to sensibly present information 
about those solutions. 

We first applied Newton's method to solve semilinear elliptic boundary value problems in |20| . 
where we sought solutions as critical points of an appropriate action functional on a suitable function 
space. This article combines the new ideas introduced in |17| . concerning nonlinear PdE on graphs, 
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with the recent advances concerning symmetry, bifurcation, and automation presented in our paper 
|19j . By automating the symmetry analysis and corresponding isotypic decompositions found in [19] 
for arbitrary graphs, we are able to apply two modified implementations of the gradient Newton- 
Galerkin algorithm (GNGA, see [20]) in order to completely automate bifurcation branch following. 
We are able to handle most difficulties that arise when encountering accidental degeneracies (see 
Definition I5.ip and high-dimensional critical eigenspaces. 

In order to catalog experimentally found solutions according to symmetry and symmetry type 
and understand the type of bifurcation that lead to the successful computation of each solution, 
we make use of the automatic generation of the information described in the bifurcation digraph 
(see [19J) corresponding to that experiment's underlying symmetry group. A visual display of 
the bifurcation digraph can also be automatically generated for human use in understanding the 
underlying variational structure and expected proliferation of solutions. 

In our bifurcation diagrams, we show plots of versus s for solutions u to Equation (TjQ) 

with parameter s. These diagrams can indicate by color or line type the symmetry of solutions, 
which is invariant on each branch, or the Morse Index (MI) of solutions, which typically changes 
at bifurcation and turning points. 

In a generalization of the notion of a contour plot, we have developed several visual represen- 
tations of solutions to the discrete nonlinear problem (fT]). To a high degree, these plots too are 
automatically generated, chosen where possible to make the symmetry of solutions visible and to 
yield a graphic that is informative and pleasing. 

Our study of the finite dimensional semilinear elliptic PdE ([I]) closely follows the related works 
concerning the PDE 

/ A« + /(k) = in ft 
[ ) \ =0 on aft, 

as well as the similar zero-Dirichlet problem; see [U1|TB] and references therein. The graph Laplacian 
L corresponds to the negative Laplacian —A from PDE theory. Both L and — A have non-negative 
eigenvalues, and both have zero eigenvalues corresponding to constant eigenvectors and eigenfunc- 
tions, respectively. Note the sign difference between Equations ([I]) and ([2]). 

Much is known about the spectrum of the graph Laplacian. See, for example, [3 t 4, 8j. Most of 
the PdE literature concerns linear problems and/or positive solutions, whereas we are interested 
in the existence and symmetry of all solutions, in particular sign-changing ones, to nonlinear PdE. 
Our first paper in this subject area [17] contains a fairly thorough list of citations relevant to the 
study of solutions to linear and nonlinear PdE, e.g., works by A. Ashyralyev, S. S. Cheng, P. G. 
Kevrekidis et al, S. T. Liu, M. Lapidus, G. I. Marchuk, C. V. Pao, Yu. V. Pokornyi, V. L. Pryadiev, 
P. Sobolevskii, J. C. Strikwerda, and G. Zhang. Applications of limiting cases where we increase 
the number of vertices and use a scaling factor to approximate solutions to nonlinear PDE on 
fractals closely follow the linear results of R. S. Strichartz and A. Teplyaev. Our survey article |16] 
summarizes some of our most relevant PDE results and provides a list of open problems in that area. 
While PdE are generated whenever PDE are discretized via finite differences on a grid, the PdE 
we study in this paper have few vertices and do not approximate PDE. Our numerical techniques, 
symmetry analysis, and existence theorems [Q[T7], apply to both PdE and PDE. By focusing on 
PdE, we can study large symmetry groups which would only arise from PDE on domains with 
dimension 3 or greater. 

In [20], the Gradient Newton Galerkin Algorithm (GNGA) was developed to investigate PDE ([2]) 
using a basis of eigenfunctions of the corresponding (continuous) linear problem to span a suitably 
large finite dimensional subspace. In [19] we adapted this algorithm in order to find many solutions 
of a PDE on a region with fractal boundary, while in [T7] we used the entire basis since n was small. 
In the current work, we modify the GNGA slightly in two different ways; the cylinder-augmented 
GNGA (cGNGA) is used to find initial solution points on new branches near bifurcation points, 
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and the tangent-augmented GNGA (tGNGA) is used to more effectively follow solution branches. 
A closely related approach to finding solutions to PDE with symmetry is found in |22l [23] . 

The paper is organized as follows. In Section [2] we handle the preliminaries, stating definitions, 
theory, and notation for graphs, graph Laplacians, symmetry, the variational method, and isotypic 
decompositions. Section [3] enumerates the various tasks we do prior to invoking the continuation 
solver which implements the GNGA to find solutions to Equation (TjQ). In particular, we discuss 
graph creation and layout from an edgelist, the computation of Aut(G), bifurcation digraphs, the 
orthonormal basis of eigenfunctions of L, and isotypic decomposition of symmetry- invariant fixed 
point spaces for bifurcation analysis. We provide some details in Section H] concerning the im- 
plementations of the tGNGA, secant method, and cGNGA. These three algorithms are used for 
following branches, finding bifurcation points, and finding new solutions on bifurcating branches, 
respectively. Section [5] outlines our algorithms and heuristics for controlling the repeated appli- 
cation of the Newton and secant code to find representative branches from every conjugacy class 
of branches. Some postprocessing details for generating contour plots and bifurcation diagrams 
are given in Section [6l Our main examples and numerical results are found in Section [7J The 
concluding Section [8] contains observations and ideas for future refinements and applications of our 
methods. 



In this section we review background and notation for graph theory, symmetry, and the GNGA. 

2.1. Graphs. Let G = (Vg,Eq) be a simple connected graph with vertex set Vq = {v\, . . . ,v n } 
and edge set Eg- The degree of a vertex Vi is denoted by d(vi). An automorphism of G is a bijection 
ol '■ Vg — > Vq such that {a(vi), a(vj)} S Eq if and only if {v j, Vj} E Eq. The symmetry group of G 
is the group Aut(G) of automorphisms. If ir is a permutation in S n then we define a n : Vq — > Vg 
by a n (vi) = iWj)- Not every permutation defines an automorphism of G but every automorphism 
of G is determined uniquely by a permutation. 

2.2. Cayley graphs. We are often interested in graphs with prescribed symmetry. We construct 
these graphs as decorated Cayley graphs [24] ■ Given a group T and a set of generators A, the 
Cayley color digraph Cay^r is a directed labeled graph (G, c) with vertex set Vg = T and edge set 



Edge (g, gd) is labeled with the color c(g, gd) = d. The group G acts on Cay A T by left multiplication; 
in fact, Aut(Cay A T) = T. To create a simple undirected graph whose automorphism group is T, 
we replace the directed colored edges of the Cayley color graph with undirected decorated edges. 
The decoration adds extra vertices along the edges. The resulting simple graph is called a decorated 
Cayley graph. 

2.3. Graph Laplacian. The Laplacian of G is determined by the matrix L defined be letting 
Ljj = d(v{), Lij = — 1 if {vi,Vj} £ Eg, and L^j = if i ^ j but {vi,Vj} Eq- This Laplacian 
can be viewed as enforcing the zero Neumann boundary condition [5]. For example, if we solve 
the appropriately scaled Equation ([T|) on the path P n for large n, we get approximate solutions 
to Equation ([2]). Consideration of other boundary conditions is possible and interesting, but is 
the subject of other and future reports. The incidence (first difference) matrix D of an arbitrary 
orientation of G satisfies L = D T D; we do not use this fact but observe that the variational 
equations for PdE most closely resemble those for PDE when expressions like Lu ■ v are replaced 
with Du ■ Dv (see [17]). The eigenvalues and corresponding eigenvectors of L are denoted by 
= Ai < A2 < • • • < A n and {ipj}j" = i, respectively. 

Let X = {(u, s) e 1" x 1 -Lu + f s (u) = 0} be the solution set of PdE Q. We write 
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where ^> m = {ipi, . . . ,ip m } is an orthonormal set of eigenvectors of L, and use the notation 
([u]* TO ,s) := (a, s) £ R m x R. Since we are working with modest sized graphs, in this paper 
we take m to be n. 

Where possible, when multiple eigenvalues are encountered, choices of the associated eigenvectors 
in ty m are made to respect symmetry in a similar fashion as was done in |18j . For details see 
Section 13.51 

2.4. Symmetry of functions. To study the symmetry of solutions to Equation ([1]) we consider 
Tq = Aut(G) x Z2, where Z2 = {1, —1} is written multiplicatively. The natural action of Tq on R™ 
is defined by 

(3) (7-u) i = /3« w -i (i ), 

where 7 = (a n ,f3) £ Tq and u £ R n . We usually write a for (a, 1) and —a for (a,— 1). The 
symmetry of it is the isotropy subgroup Sym(n) := Stab (it, To) = {7 G To | 7 ■ k = it}. Two 
subgroups Ti and Tj of Tq are called conjugate if Tj = 7r ? '7~ 1 for some 7 S IV The symmetry type 
of it is the conjugacy class [Sym(it)] of the symmetry of it. We use the notation Q := {Tq, . . . , T q } 
for the set of symmetries and S := {Sq = [Tq], . . . , S r } for the set of symmetry types. 

In general it is difficult to compute Q, but the following definition helps for some graphs. A 
generic vertex of a graph is a vertex v such that {a £ Aut(G) | a(v) = v} contains only the 
identity map. The Aut(G) orbit of a generic vertex has the same size as Aut(G). The proof of the 
following proposition follows |19j . 

Proposition 2.1. If the graph G has a generic vertex, then Q = {T < Tq \ T = Tq or — 1 T}. 

Proof. Assume that G has a generic vertex, which we label as v\. Consider the function it such 
that 111 = 1 an d itj = for i > 1. Then for any subgroup r < Tq the function 

^7-u 

has symmetry r if — 1 T, and symmetry Tq otherwise. On the other hand, only it = satisfies 
—11 = it. So, if r < Tq is an isotropy subgroup containing —1 then T = Tq. □ 

Remark 2.2. We have a counterexample which shows that the converse of Proposition 12. II is false. 
Our counterexample G is the union of a decorated Cayley graph of Z3 and a decorated Cayley 
graph of Z5, with 15 additional edges joining each element in Z3 with each element in Z5. The 
symmetry group of G is Aut(G) = Z3 x Z5. The set of symmetries Q of G consists of Tq and all 
the subgroups of Z3 x Z5, but G has no generic vertex. 

A decorated Cayley graph of any group automatically has a generic vertex, namely any of the 
vertices corresponding to elements of the group. Graphs with generic vertices are good models 
for PDE where the domain f2 has a particular symmetry. If we can find a graph G such that 
Aut(G) = Aut(O) and G has a generic vertex, then Q for G is the same as the set of possible 
symmetries of solutions to PDE ([2]) on 

We define a branch of solutions to be a maximal subset of the solution space X that is a C 1 
manifold with constant symmetry. The trivial branch {(0,s) | s £ M} contains the trivial solution 
it = 0, which has symmetry Tq if f s is odd, and symmetry Aut(G) otherwise. The positive constant 
branch is {((c, . . . , c), s) | f s (c) = 0, c > 0, s £ M.}. The negative constant branch is similarly defined. 
A bifurcation point is a solution in the closure of at least two different solution branches. We call 
the branch containing the bifurcation point the mother, and the other branches the daughters. For 
example, the (positive and negative) constant branches are daughters of the trivial branch, which 
contains the bifurcation point (0, 0) £ M. n x R. 

Equation ([1]) can be interpreted as VJ s (u) = 0, where VJ S : R n — > R n is defined by —VJ s (u) = 
—Lu+f s (u). The operator VJ S is Aut(G)-equivariant, i.e., VJ s (cm) = a-V J s (u) for all a £ Aut(G). 
Furthermore, if f s is odd, then V J s is To-equivariant. If it is a solution to Equation ([!]) with f s odd, 
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then j-u is also a solution to Equation ([I]) for all 7 G Tg. Following the standard treatment |12tll9|. 
for each Tj < Tq we define the fixed point subspace of the Tq action on V = R n to be 

Fix(ri, V) = {u G R n I 7 • u = u for all 7 G TJ. 

If / s is odd, these fixed point subspaces are VJ s -invariant. Otherwise, Fix(r"j,y) is VJ s -invariant 
for all Ti < Aut(G). Recall that a subspace W C R™ is VJ s -invariant if VJ S (W) C FY. We say 
that a subspace .4 C R n is an anomalous invariant subspace (AIS) if it is VJ^-invariant but is 
not a fixed point subspace. If A is an AIS, we sometimes say u £ A is anomalous. Note that 
the constant subspace W c := {(c, . . . , c) | c € R} is the fixed point subspace Fix(AutG, R n ) if G 
is vertex transitive. Otherwise, the constant subspace is an AIS, which we denote by A c . The 
book jllj is a good reference on invariant spaces of nonlinear operators, although our definition of 
anomalous invariant subspaces appears to be new. 

We define an anomaly-breaking bifurcation to be a bifurcation where the daughters have the same 
symmetry as the mother, and the mother is in an AIS that does not contain at least one of the 
daughters. We have not been able to describe a general theory for anomaly-breaking bifurcations. 



2.5. Isotypic Decomposition. To analyze the bifurcations of a branch of solutions with symme- 
try Ti, we need to understand the isotypic decomposition of the action of Tj on R n . 

Suppose a finite group T acts on V = R n according to the representation g i->- a g : F — > Aut(V) = 
GL n (R). In our applications we choose r G Q and the group action is the one in Equation ([3]). 
Let {a^ : F — > GL d (fc)(R) | k G K?} be the set of irreducible representations of F over M. We 

write and K when the subscript F is understood. It is a standard result of representation 
theory that there is an orthonormal basis By = Ufcex ^y ^ or ^ such that B r k ^ = IJ^j B v k ^ and 
[a g \ v {k,i)] R (fe,i) = a^ k \g) for all g G F, where vS k ' 1 ' := span(f?p fc '^). Each vS k is an irreducible 

subspace of V . Note that B^ might be empty for some k, corresponding to = {0}. The 
isotypic decomposition of V under the action of F is 

keK 

where = ©z=l Vp*'^ are the isotypic components. 

The isotypic decomposition of V under the action of each Tj is required by our algorithm. The 
decomposition under the action of Aut(G) is the same as the decomposition under the action of 
Fq. While there are twice as many irreducible representations of Fq = Aut(G) x Z2 as there are 

of Aut(G), if ofp (— 1) = I then = {0}. The other half of the irreducible representations 
have ct^\— 1) = —I. The irreducible representations of Fq and of Aut(G) can be labeled so that 

vff = V^ t{G) foTkeK Aut{G) . 

The isotypic components are uniquely determined, but the decomposition into irreducible spaces 

(k) (k) (k) 

is not. Our goal is to find B r for all k by finding the projection P r : V — > V r . To do this, 
we first need to introduce representations over the complex numbers C for two reasons. First, 
irreducible representations over C are better understood than those over R. Second, our GAP 
program uses the field C since irreducible representations over R are not readily obtainable by 
GAP. 

There is a natural action of F on W := C n given by the representation g 1— >• j3 g : F — )■ Aut(W) such 

✓tn (k) 

that j3g and ct g have the same matrix representation. The isotypic decomposition W = £B fcg ^> Wf' 
is defined as above using the set {j3^ : F — > GL-( fe )(C) | k G Ky} of irreducible representations of 
r over C. 
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The characters of the irreducible representation are X^Hd) := Trf3^ k \g). The projection 



QP : W -> W^ k) is known to be 



(4) Q? = -^rY.* (k) W 9 - 

1 1 <?er 

(k) (k) 

We are going to get the Pp s in terms of Q r s with the help of the Frobenius-Schur indicator 

" {k) :=4EX ( *V)€{-1,0,1}. 



IGI 

1 1 see 



Recall 191 that 



(i) = 1 implies = X^ k \ m which case we say is a real irreducible representation; 

(ii) z^( fc ) = implies x 7^ X^ k \ m which case we say (3^ is complex; 

(iii) z^ fc ) = — 1 implies x = X^ k \ m which case we say is quaternionic. 

Sometimes the term quasi-real is used in place of quaternionic. For k, k' £ Kp, we say k ~ k' if 
•j^w — T^( fc/ ). Complex representations come in complex conjugate pairs, so ~ is an equivalence 
relation. Then K-p can be chosen to be any complete set of representatives of the quotient set 
KyI ~. We calculate the projection operators in R™ from the projection operators in C n using the 
following formulas: 

(i) if z/ fe ) = 1 then Pf ] = Q {k) \ v and 4° = 4^ 

(ii) if = then P r fc) = {qP + Q (k) ^j \ v and Sp = 2d (k) ; 

(iii) if */(*) = -1 then P^ k) = QP \ v and 4* } = 2df\ 
for all k £ K T . 

2.6. GNGA. We now review the GNGA for PdE [IT]. Let F s : R -> R be the primitive defined 
by F s (t) = f* / s (r) dr, e.g., F s (t) = |st 2 + \t^. The actton functional J s : R n -> R is defined by 



n 



J s («) = \Lu ■ u -y]F s (ui). 

i=i 

For u,v £ R n it is easy to see that 

j;( I1 )W = -(-k + / s (tt)). ) ), 

so that n is a critical point of J s if and only if (u, s) 6 X, i.e., u is a solution to Equation (fTJ) for 
parameter s. 

For u 6 C/ m , we compute the the coefficients of the gradient vector g s (u) £ R m by 

rn 

(5) 3s(u) 3 - = Lu ■ ipj - f s (u) ■ ipj = {Ly^ a k ipk) ■ ipj - fs(u) ■ ipj = ajXj - f s (u) ■ ipj. 

k=l 

If m = n, then g s (u) = if and only if VJ s (w) = 0. When m < n, the solutions to g s (u) = are 
approximate solutions to Equation ([1]). In this paper we assume m = n, but the formulas use m 
where appropriate, keeping in mind applications to large graphs, e.g., those arising from a PDE. 
Similarly, the Hessian matrix h s (u) = (J' s '(u)(ipj, ^fc))Jj. =1 can be computed as 

(6) h s (u) jk = Ltpj • ip k - di&g(f' s (u))ipj • ip k = XjSjk - diag(f' s (u))ipj • ip k , 

where 5j k is the Kronecker delta and diag(/^(n)) is a diagonal matrix. Using the coefficient vector 
a and eigenvalues {Xj}™^ to compute the difference terms in Lu ■ ipj and Lipj ■ ip k significantly 
reduces the number of matrix and vector operations required to define the linear system for a search 
direction x satisfying h s (u)x = g s (u)- Applying Newton's method to find zeroes of (u,s) h-» g s (u) 
is the basis of our gradient Newton-Galerkin algorithms (GNGA). 
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We define the signature sig(u,s) to be the number of negative eigenvalues of the matrix h s (u) 
representing the self-adjoint bilinear operator D 2 J s (u). If (u, s) is a nondegenerate solution to 
Equation (pQ), then sig(ii, s) equals the Morse index MI(u, s). The MI can be thought of as the 
number of "down" directions of the critical point, that is, Ml(u, s) = for minima of J s , MI(u, s) = 
n for maxima, and MI(u, s) G {1, . . . , n— 1} for saddle points in between. The search direction \ can 
be found using any number of linear solvers; we use a minimum norm least squares solver to avoid 
problems with noninvertible Hessians h s (u). Noninvertible Hessians inevitably occur at bifurcation 
points, and fold points (points where the solution branch is not monotonic in s). When the Hessian 
is singular, the eigenspace of the Hessian with eigenvalue is called the critical eigenspace, and is 
denoted by E. 



In this section we describe the various tasks that must be performed prior to approximating 
solutions to Equation (UJ. In particular, we must create the edgelist, visualize the graph, analyze 
the symmetry of the problem, compute the possible bifurcations, and generate suitably bases from 
the eigenfunctions of the Laplacian. The data files generated during preprocessing are used by 
the continuation solver, as well as in the postprocessing phase when creating graphics in order to 
visualize the results. All of these files for a single graph, Example 17. 2} can be found at the website 

\protect\vrule widthOpt\protect\href {http : //NAU . edu/ Jim. Swif t/PdE>{http : //NAU . edu/Jim . Swif t/P< 

3.1. Graph Creation. A graph G is determined by an edgelist file. Each line contains a pair of 
integers i and j, indicating that {v j, Vj} S Eq. This file is usually created by a text editor. We also 
have the option to create the edgelist file automatically by GAP [13] as a decorated Cayley graph 
of a given group (see Section [7. 6p . This is the main human input for our process. 

3.2. Graph Layout Code. To create a visualization of the graph, we use a standard spring 
embedding algorithm to create an embedding £ : Vq — > M? of the graph. This is done by a C++ 
program. The program starts with a random placement of the vertices, that is, I is initialized with 
random values. We then calculate the "force" Fi = Ei + Hi on each vertex Vi, where E{ is generated 
by an equal "electric charge" Q on the vertices and Hi is generated by "springs" of natural length 
v replacing the edges of the graph. Specifically, with dij = i{vi) — i(vj), we have 



Experiments show that the values Q 2 = 1, e = 0.001, v = 1 and D = 1.1 work well. Iteratively 
replacing £{vi) by £(vi) + SFi using a stepsize of 5 = 0.1, we simulate the movement of this physical 
system with added damping until an equilibrium is reached. This stable position usually shows some 
aspects of the symmetry of the graph. The complexity of the layout is defined to be the number 
of distinct distances between vertices. The program tries several initial positions and picks the 
layout that minimizes the complexity. Layouts with higher complexity are also stored for possible 
use. Finally, we rotate the optimal placement so that the most common edge slope is horizontal. 
The output of the program is a file containing the coordinates of the vertices. We create figures 
automatically from this file using Gnuplot, )$^-pic and Mathematica, together with solution data 
generated by the continuation solver. 

3.3. Automorphism Group Code. To analyze the symmetry of the graph we need to find its 
automorphism group Aut(G). This is done by Nauty |15| . which is a very efficient program that can 
handle fairly large graphs. It creates a file containing all the permutations or only the generators 
of the automorphism group. This file is the input of the GAP program that performs the full 
symmetry analysis. 
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3.4. Symmetry Analysis Code. In this subsection we give some details of the GAP computations 
done to analyze the symmetry of the problem. Since GAP uses irreducible representations over the 
complex numbers, some care must be taken. 

To compute the set of symmetries Q, we need the following definition. If T acts on V and U is 
a subspace of V then 

pStab([/, r) = {7 e T | 7 • ii = u for all u G U}. 

The isotropy subgroups of the Tq action on R n are precisely the subgroups T of the finite group Tq 
which satisfy 

P stab(Fix(r,v),r) = r. 

This computation is easily performed by GAP. 

Next, we determine the possible symmetries of the daughters of a bifurcation point with symme- 
try G Q. For each i, we use GAP to find the irreducible representations {a^ \ k E -Kr\} of Tj, 
and the characters ■ The characters are used to produce the projection operators Qj?) defined 

/TTh • (k) 

in Equation ([4]). The isotypic components Vp. of the Tj action on W = C n are computed as the 
row spaces of the projection operators. The kernel of the irreducible representation, denoted T' ik , 

(k) 

is also computed by GAP. Then for each k 6 for which Vp. is nontrivial we generate the set 

(k) 

Hi^k of isotropy subgroups of the Tj action on Vp. . The set H^k is partially ordered with Tj at 
the top and T' ik at the bottom, and is often called the lattice of isotropy subgroups [T2J [19]. If 
there are no subgroups in %i : k properly between T{ and Tj for some Tj £ Hik then Tj is called 
a maximal isotropy subgroup. For each of these maximal isotropy subgroups there is a possible 
generic bifurcation from a mother with symmetry Tj to a daughter with symmetry Tj, represented 
by the bifurcation arrow 

The bifurcation arrows always join isotropy subgroups in Q, since H i k C Q. For each of these 
bifurcation arrows, GAP computes the groups Ti/T^j, Nr^Tj), and Nr i (Tj)/Tj (the normalizer of 
Tj in Ti, denoted Nr^Tj), is the largest subgroup of Tj for which Tj is a normal subgroup). 

At a nondegenerate bifurcation of a solution with symmetry Tj, the critical eigenspace E is an 

(k) 

irreducible subspace lying in one of the Vp and we say that the mother undergoes a bifurcation 

with Ti /T'i k symmetry. Note that Tj /T'- k acts freely on E C Vpv ■ The bifurcation arrows with a 
given label k represent the symmetries Tj of the daughters that are expected to bifurcate from the 
mother when E C V^ k \ If the system ([1]) is restricted to Fix(Fj,]R n ) then the effective symmetry 
of the bifurcation is N^.(Tj)/Tj. For example, if N^.(Tj)/Tj = Z2, then there is a pitchfork 
bifurcation creating two conjugate daughter branches. For details, see [T2l [19] . 

We say that two arrows Tj — — > Tj 1 and Tj — — > Tj 2 are equivalent HTj 1 and Tj 2 are conjugate. 
Since we only seek non-conjugate solutions in our continuation solver, one bifurcation arrow from 
each equivalence class, together with its auxiliary information, is written into a file by GAP. This 
file is used by our continuation solver when a bifurcation is encountered, as described below. 

The amount of material contained in the bifurcation arrows is overwhelming. To summarize the 
collected information about the possible bifurcations we draw a bifurcation digraph (see [19], which 
gives an equivalent definition). This is an extension of the usual lattice of isotropy subgroups. For 
an example, see Figure[2j Our continuation solver requires the label of the irreducible representation 
k, but does not require information about the group Ti/T' ri k ; on the other hand, humans find the 
symmetry group of the bifurcation more informative than k, so it is included in the bifurcation 
digraph instead of k. 

Definition 3.1. The bifurcation digraph of the Tq action on a real vector space V = W 1 is a directed 
graph with labeled arrows between the symmetry types. We draw an arrow from \Tj\ to [Tj] if and 
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only if r,- is conjugate to a maximal isotropy subgroup of the Tj action on some isotypic component 

■ The label on this arrow is Tj /T^ k , where 1^ k is the kernel of the Tj action on Vp . We use 
the arrow types 

solid ► if N Ti (T j )/T j ^Z 2 , 

dashed - - -» if Nr^T^/Fj ^ Zi, and 
dotted > otherwise, 

to indicate the nature of the bifurcation. 

(k) 

Assume that E C Vj-i. is an irreducible subspace, and a critical eigenspace of a bifurcation 
point with symmetry Tj. The group N^^T^/Tj acts freely on Fix(F j,E). The size of this factor 
group, which determines the arrow type in the bifurcation digraph, is passed to the continuation 
solver. The arrow type gives us information about the dimension of Fix(Tj,E). For real vector 
spaces V = M n , the solid and dashed arrows imply that dim^ Fix(Fj, E) = 1, whereas the dotted 
arrows give dim^ Fix(Tj, E) > 1. The first two arrow types, with 1-dimensional fixed point spaces, 
are called EBL bifurcations since the Equivariant Branching Lemma [12] guarantees (under certain 
conditions) that solutions with symmetry Tj are born at the bifurcation. These bifurcating branches 
are particularly easy to follow numerically, since there is only one critical eigenvector (up to a scalar 
multiple) with the symmetry Fj. The EBL can be considered an extension of the classic bifurcation 
results from |21j . 

If there is a dotted arrow to Tj and certain nondegeneracy conditions hold, then there is a 
daughter with symmetry Tj born for gradient systems |12| . No general theory predicts where the 
daughters lie when projected to Fix(F j,E); our approach to following such branches numerically 
requires randomly choosing perturbations within Fix(Fj,E). 

The bifurcation digraph is often very complicated so we use condensation classes instead of con- 
jugacy classes to get a simpler condensed bifurcation digraph. See Figure [2] for an example. Let 
(p G Aut(ro). If 4>(Fi) is a symmetry group for all symmetry groups Tj then we say that <j) is symme- 
try preserving. The symmetry preserving automorphisms form a subgroup Aut c (ro) of Aut(Fo). A 
symmetry preserving map induces a permutation of symmetries so Aut c (ro) acts on Q. Conjugate 
elements of Q are on the same orbit since the inner automorphisms of Tq are in Aut c (ro). Hence 
Aut c (ro) also acts on S. The orbit equivalence classes of this latter action are called condensation 
classes, and are computed automatically by a GAP program. The condensed bifurcation digraph 
is the quotient of the bifurcation digraph by the orbit equivalence of the Aut c (ro) action on S. 
Hence, the vertices of the condensed bifurcation digraph are the condensation classes. 

3.4.1. Digraph Layout Code. Using a C++ program similar to the graph layout code, we generate 
a layout for the bifurcation digraph and the condensed bifurcation digraph The difference is that 
the vertices can move horizontally but their vertical position is determined by the size of the group 
they represent. Graphics of the layouts are then created by Gnuplot and }^-pic. 

3.5. Basis Generation Code. To generate ^> m = {ipi, ■ ■ ■ ,ip m }, we use a somewhat complicated 

k 

procedure that increases the efficiency of the GNGA. We first pick a single arrow Tq ► 

for each k G Ky q . Using Mathematica, we then find a basis of eigenvectors of the symmetric 
operator L restricted to the invariant subspace nV^. We compute the orthonormal basis 

of Vp ( o fc) using the Gram-Schmidt process on {7 • v \ 7 G Fq,v G D^}. We start the Gram-Schmidt 

process with the already orthonormal set JDW so that these elements survive in the resulting basis. 
The eigenvalues of L and the corresponding eigenvectors in *$> m = U^^^^ are written to a file. 

(k) 

The basis generation code also calculates the projection operator Pp. for each i and k G K?. 
from the characters and permutations produced by our automorphism group and symmetry analysis 
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(k) i /• i 

codes. The isotypic component Vp. is the range of -PA • For each i and k the coordinates in \E' 



(k) 



of the elements of are written to a file. 



I',; 



4. Newton's Method with Constraints 

To follow branches or find new bifurcating branches, we treat the parameter s as the (m + l) st 
unknown. Thus, when we say p = (a,s) G M m+1 is a solution, we mean that u = ^2ajipj solves 
Equation ([I]) with parameter s. We restrict the search for a particular solution to some hypersurface 
in R m+1 , satisfying an (m+l) st equation of the form /t(a, s) = 0. In this section we will describe two 
different choices of At, one for following branches and another for finding new branches emanating 
from bifurcation points. In the first case, we take an old and current pair of solutions p Q and p c 
along a symmetry invariant branch to obtain a reasonable initial guess p g for iterating to find a new 
solution p n satisfying the constraint of lying on a hyperplane normal to the branch. The constraint 
used at a bifurcation point p* instead forces the new solution p n to have a projection onto the 
critical eigenspace of a specified norm. 

In either case, the iteration we use is: 

• compute the constraint k, gradient vector g := g s (u), and Hessian matrix h := h s (u) 

h % 

(Va*) T fe 

-(a,s)-x, u = Y, a j^j- 
Equations ([5|) and ([6]) are used to compute g and h. The (m + l) st row of the matrix is defined 
by (V a K, |j) = Vk G M m+1 ; the search direction is % = (Xa,Xs) G M m+1 . Since this is Newton's 
method on (g, k) G R m+1 instead of just g G R m , when the process converges we have not only that 
g = (hence p = (a, s) is a solution to Equation ([I])), but also that k = 0. 



solve 



(a,s 





Xa 
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. Xs _ 




K 



4.1. Tangent-augmented Newton's method (tGNGA). We use the tGNGA to follow 
branches. Given two consecutive solutions p Q and p c on a given branch, we compute the (ap- 
proximate) tangent vector v = (p c — p )/\\Pc ~ Po\\ G M m+1 . The initial guess is then p g =p c + cv. 
In our experiments the speed c has a minimum and maximum range, for example from 0.01 to 0.4, 
and is modified dynamically according to various heuristics (see for example Figure H]) . For the 
tGNGA, the constraint is that each iterate p = (a, s) must lie on the hyperplane passing through 
the initial guess p g , perpendicular to v. That is, n(a,s) := (p — p g ) ■ v. Easily, one sees that 
(V a K(a, s), §f (a, s)) = v. In general, if f s has the form f s (u) = su + H(u), then || = —a. For 
example, when f s {u) = su + n 3 , a calculation shows that g s (a)j = o,j(Xj — s) — (J2T=i a ktpk) 3 ■ ipj, 
hence ^| = —a. Newton's method is invariant in this plane so that in fact x • ^ = at each step. 
Hence, the linear system to be solved each iteration can be described by: 



h |2 " 

as 




Xa 
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_ {v a ) T Vs 




. * s _ 








where v = (v a ,v s ) G R m+1 . 
satisfying the constraint. 



Our function tGNGA(p 9 ,w) returns, if successful, a new solution p n 



4.2. The secant method and processing bifurcation points. In brief, when using the tGNGA 
to follow a solution branch and the MI changes at consecutively found solutions, say from k at the 
solution p to k + d at the solution p c , we know by the continuity of D 2 J S that there exists a third, 
nearby solution p* where h is not invertible and the r th eigenvalue of h is zero, where r = k + [^] • 
Let po = p , p\ = p c , with /3o and /3i the r th eigenvalues of h at the points po and p±, respectively. 
We effectively employ the vector secant method by iterating 
(Pi -pi-\)Pi 



Pg = 
Pi+i 



Pi 



tGNGA(p g , 



v) 
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until the sequence (pj) converges. The vector v = (p c — p )/\\Pc ~ Poll is held fixed throughout, 
while the value is the newly computed r th eigenvalue of h at p^. If our function secant(p ,p c ) 
is successful, it returns a solution point p* = (a*,s*), lying between p Q and p c , where h has r 
zero eigenvalues within some tolerance. We take the critical eigenspace E to be the span of the 
corresponding eigenvectors. If p* is not a turning point, then it is a bifurcation point. 

4.3. Cylinder-augmented Newton's method (cGNGA). The cGNGA is used to find initial 
solution points on new branches near bifurcation points p*. After such a point has been detected and 
the corresponding critical eigenspace has been computed, we search for a new solution bifurcating 
from the main branch by running the cGNGA. The first step is to choose a subspace E of the 
critical eigenspace. 

To ensure that we find the mother solution rather than a daughter, we insist that the Newton 
iterates belong to the cylinder C := {(a, s) : ||Pg(a — a*)|| = e} , where Pe is the orthogonal 

projection onto E and the radius e is a small fixed parameter. The input file to the continuation 
solver sets the value of e. At a symmetry breaking bifurcation the critical eigenspace is orthogonal 
to the fixed point subspace of the mother, so the mother branch does not intersect the cylinder. 
We conjecture that at anomaly breaking bifurcations the critical eigenspace is orthogonal to the 
AIS of the mother branch. It is not true that the critical eigenspace is orthogonal to the mother in 
general, although we observed this to be true for all the numerical results we include in this paper. 
Consider Equation (fTJ) with f s {u) = (su + u 3 )(u 2 — 1), which has a bifurcation where two solution 
branches in A c cross at u = (1, 1, . . . , 1), s = —1. At this bifurcation point a* and the critical 
eigenvector are parallel. 

The constraint we use to put each Newton iterate on the cylinder is n(a, s) = ^(\\Pe{cl — a*)\\ 2 — 
e 2 ) = 0. The initial guess we use is p g := (a*, s*) +e(e, 0), where e is a randomly chosen unit vector 
in E. Clearly, p g lies on the cylinder C. A computation shows that V a «(a, s) = Pe(o- — a*), and 
§f (a, s) = 0. Hence, the search direction x is found by solving 



h |f 
{P E (a-a*)) T 





Xa 
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. Xs _ 




K 



Again, ^f = —a when f s has the form f s (u) = su + H(u). When successful, cGNGA(p*,p 9 , E) 
returns a new solution p n of Equation (fTJ) that lies on the cylinder C. 

We take E to be various low-dimensional subspaces of the critical eigenspace, corresponding to 
the possible symmetries that bifurcations theory predicts must exist. For example, at an EBL 
bifurcation E is spanned by a single eigenvector. When the dimension of E is greater than one, 
we call cGNGA repeatedly with several random choices of the critical eigenvector e. The details 
are given in Equation ([7]) and Algorithm [3l The theory we apply does not guarantee a complete 
prediction of all daughter solutions. Therefore we also call cGNGA with E equal to the full critical 
eigenspace. In this way, if the dimension of the critical eigenspace is not too big we have a high 
degree of confidence that we are capturing all relevant solutions, including those that arise due 
to accidental degeneracy and that are neither predicted nor ruled out by understood bifurcation 
theory. 

5. Continuation Solver 

Our continuation solver is implemented in C++. We start our search for solution branches with 
the known solution p c = (0, sq) G M m+1 , which lies on the trivial branch, together with the initial 
direction vector v = (0, ±1) G R m+1 , which points in the direction of another solution on the trivial 
branch. The branch queue is initialized with the job (p c , v). Every job in the branch queue is fed to 
f ollow_branch until the branch exits some window in M m+1 . After every new point is computed, 
f incLbif points is called. If a bifurcation point is found, f incLdaughters is called and a job is 
added to the branch queue for every solution found. For efficiency, f incLdaughters only returns 
solutions on distinct, non-conjugate bifurcating branches. The process stops when the branch queue 
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while (p c G window and c > r) 

set p g = p c + cv 
set p n = tGNGA(p 9 , u) 
if p n is unacceptable 
halve speed c 



else 



use heuristics to adjust speed c 
set pc = p c 
set p c = p n 



-Po 



set v = m M 

||pc-Po|| 

forall p* G f ind_bif points(p G , p c ) 

compute critical eigenspace Ei of p* 
forall q{ G f ind_daughters(p|, E{) 

* J 

l|p*-9ill 

add (p*,vj) to branch_queue 



Algorithm 1: f ollow_branch(p c , v) 



if |MI( Po )-MI(p c )| = 

return { } 
set p = secant(p D ,p c ) 
compute critical eigenspace E of p 
if dim(£) = | MI(p ) - MI(p c )| 

if p is not a fold point 
return {p} 

else 

set p g = (p + p c )/2 



set p = tGNGA(p g , i>) 

return f ind_bifpoints(p Q ,p) U f ind_bif points(p,p c ) 
Algorithm 2: f ind_bifpoints(p ,p c ) 



is empty. Thus, our continuation solver finds a representative branch from each conjugacy class of 
branches connected to the trivial branch, within the chosen window. 

5.1. Branch Following. Once an instance of f ollow_branch has started, consecutive solutions p Q 
and p c are used to generate the next direction vector v = i i^ c ~^° ii . The speed c is modified accord- 
ing to the scale and complexity of features, e.g., severe turning points, proliferation of proximal 
bifurcation points, or failure of the algorithm to converge. When the algorithm converges especially 
quickly, for example, we use other heuristics to increase the speed as our guesses are somehow too 
good. In that way, many solution points are found near trouble spots while much fewer are needed 
on long, featureless parts of a branch (see for example Figure Q]). Each point along the branch 
together with its corresponding data is written to a file, to be used later in generating bifurcation 
diagrams and reports, as well as for diagnosing the occasional failure. Algorithm Q] is executed 
repeatedly until the branch queue is empty. 

The implementations of the tGNGA and secant methods are straightforward following Section [5J 
Details concerning find_bif points and f ind_daughters can be found in Algorithms [2] and O 
respectively. 
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S 



FIGURE 1. Bifurcation diagram for the first primary branch of the Cayley graph 
of S3 (see Section I7.4H . The graphic demonstrates how the density of computed 
points (denoted by crosses) is increased near interesting features. We use heuristics 
to adjust the speed. For example, the speed c is halved if tGNGA fails to converge 
in four iterations. Further, the speed is multiplied by a factor in (0, 2] based on the 
angle formed by the last three points, where the factor is 1 if the angle is 0.1 radians. 



5.2. Finding bifurcation points. When a MI change is observed between two consecutive solu- 
tions p and p c on a branch, we know that there exists at least one degenerate intermediate solution 
p* on the branch. Since there can be multiple degenerate points on the branch in the branch seg- 
ment, our bifurcation point finding algorithm recursively calls the secant method of Section T4.2I (see 
Algorithm [2]) . 

5.3. Finding daughter branches. Given a bifurcation point p* = (u*,s*) with symmetry Tj 
and its critical eigenspace E, we want to find all of the bifurcating branches. Our function 
f incLdaughters, described in Algorithm [3l finds as many daughters as it can. 

If p* is nondegenerate, then E is an irreducible subspace so it is contained in exactly one iso- 

(k) . (k) 

typic component Vp\ of the Tj action on V = M. n and has dimension d r . . The first step of 
f incLdaughters is to find the intersection of E with each isotypic component. These determine 
the set K = {k | E n V^ k) ^ {0}}. Here we use the set of bases {B^ \ k 6 -Ki\} computed in 
Section 13.51 

A bifurcation point is degenerate if E is contained in the fixed point subspace Fix(rj,V) of 
the mother. When u* = 0, E cannot be contained in the zero-dimensional fixed point subspace 
Fix(ro, V) = {0} of the mother. In all other cases, we label the irreducible representations of Tj so 

that the trivial representation is k = 1, and Fix(Fj, V) = Vp . 
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compute K and J for Ti and E 
compute Ej for j G J 
set list_of .daughters = { } 
if dim(£) > 1 
set E = E 
set J = J U {0} 
forall j € J 

set num_no_changes = / nc (dim(.Ej)) 

set no_changes = 

while no_changes < num_no_changes 

choose a random e G -Ej with ||e|| = £ 

set p g = (u* + e, s*) 

set g = cGNGA(p*,p 9 , Ej) 

if the Tj orbit of q and list_of .daughters are disjoint 
set no_changes = 
add q to list_of -daughters 

set p g = (u* — e,s*) 

set q = cGNG k(p*,p g ,Ej) 

if the Ti orbit of q and list_of .daughters are disjoint 
add q to list_of -daughters 

else 

increment no_changes 
return list_of .daughters 

Algorithm 3: find_daughters(p*,.E) 

Definition 5.1. The bifurcation point p* with symmetry Ti has accidental degeneracy if any of 
the following conditions hold: 

(1) K contains 1 and Tj ^ Tq; 

(2) K is not a singleton set; 

(3) dim(£ n Vff) > 4*? for some k G K. 

Then we say that p* has a degeneracy of Type 1, 2 or 3 respectively. 

Condition (1) says that E has a nontrivial intersection with the fixed point subspace of the 
mother. If E is not an irreducible subspace then condition (2) or (3) holds. 
The set of expected bifurcating symmetry indices of the daughter solutions is 

k 

J := {j | Ti > Tj is a bifurcation arrow for some k 6 K or (j = i / and 1 G K)}. 

For each j S J, we define 

(7) E j = EnFix(T j ,V). 

For nondegenerate bifurcations, Sym(ii* + e) = Tj for all nonzero e G Ej, since is a maximal 
isotropy subgroup of the Ti action on E. Our algorithm uses initial guesses u* + e with ||e|| = e in 
the cylinder-augmented Newton's method function cGNGA to look for solutions with symmetry Tj. 

The daughters are found by repeated applications of cGNGA. A heuristic function f nc : N — > N 
with /nc(l) = 1 takes as input the dimension of Ej or E, and outputs the number of cGNGA 
consecutive calls allowed without finding a new daughter. The default function is defined by 
fnc(d) = 1 + 20(d — l) 2 , but this can be changed if one suspects that a daughter branch has not 
yet been found. 

The f ind_daughters subroutine prints information such as the number of random choices it 
took to find a new daughter, so that the user can modify the f nc function if desired. 
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The search for solutions in all of E when dim(-B) > 1 is included to find possible daughter 
branches with symmetry not predicted by any bifurcation arrows. This is needed at bifurcation 
points with Type 2 degeneracy, as in Example 17.41 Daughters with submaximal symmetry exist 
for bifurcations with certain symmetries [10] even when E is an irreducible subspace, although we 
did not encounter this in the examples we studied. 

6. Postprocessing 

Since even a small graph can have a large symmetry group and other features which lead to a pro- 
liferation of solutions via possibly complicated bifurcations, we must artfully display select subsets 
of our results in a human understandable format. This section briefly describes the methodology 
and tools we have developed which process the output from the preprocessing and continuation 
solver phases in order to generate graphics automatically, edit and annotate those graphics, and 
research new variational and symmetry phenomena. 

Our heuristics for automatically changing speed and retrying Newton's method with better initial 
guesses are sufficient to generate all the results presented in Section and many more. In a few 
instances some adjustment of the initial speed was required. This adjustment is facilitated by 
files which track every solution and present all associated information in human readable formats. 
Generally, these files also contain the actual data that the programs described in Sections 16.11 and 
16.21 use to generate graphical output. 

6.1. Contour plots. Solutions in X are displayed with a contour plot program written in Math- 
ematica. The contour plot program uses the embedding of the graph found by the layout program 
described in Section 13.21 The vertex Vi is colored white if ui > 0, gray if Ui = 0, and black if 
Ui < 0. Furthermore, the vertex is shown as a disk whose area is proportional to \v,i\. If \ui\ is 
below some cutoff, a small disk is drawn. When a solution u £ W 1 is passed to the contour plot 
program, each of the solutions in the group orbit {7 • u | 7 G Aut(G)} are tested by heuristics that 
attempt to find which one is the best. We plot the solution u that minimizes the size of the set 
{itjWj||(ijj|| I 1 < i < j < n}. To break a tie, the program chooses a solution which has a horizontal 
or vertical line of reflection symmetry, if such a solution exists. 

We say a symmetry of a solution is visible if it is also a symmetry of the contour plot. By 
changing several parameters the layout program can easily be made to generate alternate layouts 
which may make more symmetry of a given solution visible. Layouts can also be entered by hand 
or copied from the output of other programs. Once we have viewed a solution's contour plot for a 
particular layout, we can view and save any solution in the orbit of that solution. Saved graphics 
are most easily viewed in an automatically created HTML file which annotates each representative 
solution with useful information such as Morse index, symmetry and symmetry type, J value, 
branch number, and bifurcation history. In these ways, we greatly reduce the human effort needed 
to generate informative graphics in a format suitable for publication. 

6.2. Bifurcation diagrams. A (schematic) bifurcation diagram is the graph of {(s, y(u)) | (u, s) E 
X} where y : M. n — > M is some schematic function |12j . The schematic function y is needed to 
reduce the graphics to two dimensions. A good choice such as the taxicab norm y(u) = \\u\\i = 
\u\ \ + \u2\ + • • • + \u n \ visually separates branches. The To-invariance of this choice ensures that 
only one curve is shown for each equivalence class of solution branches. Thus, it avoids apparent 
discontinuities in the diagrams due to inconsistent choices of representatives of orbit classes for 
bifurcating branches. 

In |19j we used the value of a PDE solution a at a generic point of the domain. We actually 
solved a PdE, and defined y by y(u) = U{ for a fixed i, where the vertex Vi has trivial symmetry 
as described in Proposition 12.11 Note that for graphs in general (for example cycles), there may 
not be a generic vertex. While this choice of y is not a To-invariant function, we were able to get 
meaningful bifurcation diagrams without redundant branches by exploiting the simplicity of the 
symmetry group H)q in a way that cannot be done for general groups [18J. 
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In the current project, we also investigated schematic functions of the form y(u) = y w (u) = 
Yli=i w i\ u i\> f° r some choice of weight vector w. We only include results using w = (1,...,1), 
which makes y(u) = \\u\\i, but find this topic a interesting area for future research. 

7. Examples 

We considered many different graphs in our numerical experimentation, with an eye for examples 
that revealed interesting phenomena in symmetry, bifurcation, or variational structure. Among 
other things, we want to know which of the possible symmetries are represented in the solution 
space X, how the symmetry of solutions relate to the symmetries of eigenfunctions of the Laplacian, 
and what are the relationships between Morse index and nodal structure. These questions were first 
raised for PdE in |17| . We also investigate several examples of anomaly breaking bifurcations. The 
chosen examples demonstrate capabilities such as our ability to handle high multiplicity bifurcations 
and accidental degeneracies. In general, the output automatically generated during our experiments 
was sufficient for the creation of the graphics and tables included in this section. The following is 
an index of the experiments that we have decided to include in this section. 

17.11 The path P3. We demonstrate the continuation solver for five different nonlinearities, not 

all odd nor all superlinear. Our code works without modification when f s is not odd. 

We discuss the branch of constant solutions present in all our experiments. An accidental 

degeneracy of Type 1 is featured. 
17.21 The cycle C4. Branches connected to the trivial branch and the existence of solutions of 

every possible symmetry are discussed. 
17.31 Graphs with no symmetry. The smallest graphs with no symmetry have 6 vertices, and 

there are 9 such graphs. We show results for the two graphs that have AIS other than 

A c - These AIS are associated with integer eigenvalues of L and lead to anomaly-breaking 

bifurcations. 

17.41 A decorated Cayley graph of the symmetric group S3. We demonstrate that we can automat- 
ically generate a graph and solutions to Equation ([1]) on that graph with a predetermined 
symmetry group. We highlight an accidental degeneracy of Type 2 at an integer eigenvalue. 

17.51 A decorated Cayley graph of Z5. This example has several non-EBL bifurcations. The 
critical eigenspace for the bifurcations with Z10 symmetry can be in either of two different 
isotypic components; the same is true for the bifurcations with Z5 symmetry. We show 
contour plots of the bifurcating solutions that occur in the different cases. 

17.61 A decorated Cayley graph of the quaternion group Q. We construct an example with several 
occurrences of bifurcations with Q symmetry. At these non-EBL bifurcations, all points in 
the 4-dimensional critical eigenspace (except for the origin) have the same symmetry. 

17.71 The Petersen graph. Some information concerning the number of Newton iterations and 
overall computing time is presented for this fairly complicated example which has bifurcation 
points of high multiplicity. We find a numerical counterexample to a variational and nodal 
structure conjecture. 

17.81 The dodecahedron. In this example we encounter an accidental degeneracy of Type 3. The 

degeneracy is explained by an AIS that relates certain solutions on the dodecahedron to 

solutions on the Petersen graph. 
17.91 The truncated icosahedron (soccer ball). This example has more vertices and a large number 

of high multiplicity eigenvalues. We choose a layout that visually shows the resemblance of 

several solutions to spherical harmonics. 

7.1. The path P3. In this very simple example one can easily see the entire bifurcation digraph 
and possible symmetries of solutions. It is an easy exercise to increase the number of vertices in the 
path and consider scaling in order to approximate solutions to an ODE with Neumann boundary 
conditions. 
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Table 1 . Symmetries for the path P3 (Example IT. ip . The first column shows the 
isomorphism class of the elements in a condensation class; the second and third 
columns give the symmetry type and the symmetry. The fourth column shows 
contour plots for selected solutions with each symmetry type. Two solutions with 
symmetry T\ are shown, the constant solution on the left is in the AIS A c . 




5*o 



Si, S* 2 



S 3 



Figure 2. The bifurcation digraph (left) and condensed bifurcation digraph (right) 
for the Z 2 x Z 2 action on P3 (Example I7.ip . The elements in each condensation 
class are enclosed in a box. The small numerals on the arrows tell the number of 
connections emanating from each symmetry type in a box. A missing small numeral 
means 1. 



Let G = P3 be the path with three vertices. Then Aut(G) = (an 3)) = Z 2 , and so To = Z 2 x Z 2 . 
There are four symmetries in Q, shown in Table [TJ The symmetry types are all singletons, with 

s t = [FA = {r 4 }. 

The bifurcation digraph and condensed bifurcation digraph (see Section 13.4ft are shown in Fig- 
ure [2j The automorphism group Aut(ro) is isomorphic to Z 2 , and is generated by <j), where 
( K Q! (i3)) = — a (i3) an d (j>{— 1) = —1- Thus, <p interchanges Ti and T 2 , while leaving Tq and T3 
fixed. There are three condensation classes, as seen in Figure El 

Figure [3] shows the bifurcation diagrams for several nonlinearities that can be chosen by a flag in 
our continuation solver; our implementation handles non-odd nonlinearities with no modification, 
provided / s (0) =0. In addition to the functions shown in Figure El our code works for many 
other families of functions. For example, we can use f s (u) = sinh(su), which is not of the form 
f s (u) = su + H(u) since it has a nonlinear dependence on s. We also performed experiments with 
asymptotically linear nonlinearities. The solutions which bifurcate from u = 0, s = in Figure[3]are 
all constant solutions of the form u = (c, . . . , c). These constant solution branches satisfy / s (c) = 0. 
For example, in the first diagram, with f s (u) = su + u 3 , the constant solution has c = y/—s. The 
Hessian evaluated at a constant solution is a diagonal matrix with ha = \ — f' s (c) . It is an exercise 
to determine the values of s at which the Hessian is singular on this branch. 

As noted in Section \2A\ we find that W c = A c is an AIS since P3 is not vertex transitive. In the 
first diagram in Figure [3l the constant branch has a bifurcation at s = —1.5 to a daughter u A c . 
This is not a symmetry-breaking bifurcation since both mother and daughter have symmetry type 
S\ (see Table [1]). The critical eigenspace lies in the fixed-point subspace of the mother, hence this 
bifurcation point has accidental degeneracy of Type 1 (see Definition 15.1 p . 
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Figure 3. Bifurcation diagrams for the graph P3 with various nonlinearities (Ex- 
ample I7.ip . The first diagram is for our standard odd, super linear nonlinearity. 
All of the diagrams use the taxicab norm [|it||i plotted against s. Note that the 
nonlinearities featured in the bottom row are not odd, but our procedures still 
work. Extending results from |17j . one computes that the secondary bifurcations 
on the constant branch for the five cases respectively are at: s = — nonexistent; 
s = —2Xi; s = — for c < and s = — Aj for c > 0; s = — Aj for c > and 
nonexistent for c < 0. 
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Figure 4. Bifurcation diagram for C4 (Example I7.2p . showing all branches con- 
nected to the trivial branch. The symmetry type of each branch is indicated. There 
are 11 possible symmetry types for this system, but no solutions with symmetry type 
fig are found. One more branch (not connected to the trivial branch) is found in |14| . 
but this branch does not have symmetry type 5s either. The output files automat- 
ically generated by our suite of programs for this graph, including this bifurcation 
diagram, can be viewed at |http : / /NAU . edu/ Jim . Swif t/ PdE, 

7.2. The cycle C4. We investigated far too many families of graphs to include them all, but make 
a brief mention of the cycle C4 due to two interesting phenomena that we observed. Figure H] shows 
every branch that is connected to the trivial branch. Lee and Neuberger [14J found one additional 
branch for C4 that is not connected to the trivial branch, hence is missing from Figure HI The cubic 
system (JTJ) with the default nonlinearity has at most 3 4 = 81 real solutions. Lee and Neuberger 
found exactly 81 real solutions for s < s* ~ —3 by using their asymptotic form of solutions for 
large, negative s. 

Also, notice that there is no branch of symmetry type Sg in our figure. With the vertices 
numbered cyclically, functions of symmetry type Ss are of the form u = (a, b, —a, —b), with a ^ b 
nonzero real numbers. The additional branch found in p3] does not have symmetry type Ss either, 
which provides strong evidence that some systems do not realize all possible symmetry types. 

7.3. Graphs with no symmetry. We considered graphs with no symmetry. In this case, the 
set of possible symmetries of functions is Q = {ro,ri}, where T$ = Z2 and T± = {1}. The 
trivial solution has symmetry Tq and all other solutions have trivial symmetry. Thus, nontrivial 
solutions cannot undergo symmetry-breaking bifurcations. The "expected" behavior based solely 
on symmetry theory is to have Z2 bifurcations at the eigenvalues of L, which are typically simple, 
with no secondary bifurcations. As in Example 17. 1\ there are anomaly-breaking bifurcations of 
constant solutions in A c . In this section we describe another AIS that is present for some graphs 
with trivial symmetry. 
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Figure 5. Sign-changing anomalous solutions on two nonsymmetric graphs. In 
both cases the solutions lie on a primary branch bifurcating at s = 3. All secondary 
and tertiary bifurcations are anomaly-breaking for these two graphs. The bifurcation 
diagram for the graph on the left is shown. 



We have done automated experiments computing the automorphism groups of all graphs with 6 
vertices or fewer. Other than the graph with one vertex and no edges, all graphs G with 5 or fewer 
vertices had nontrivial Aut(G). There are exactly 9 graphs with 6 vertices and Aut(G) = {1}. Two 
of these graphs have AIS other than A c . Non-constant anomalous solutions to Equation ([1]) for 
these two graphs are shown in Figure [5j For both of these graphs, the vertices can be numbered so 
that the AIS is 

A 2 = {(a, a, b, b, b, b) \ a G R, b € R}. 

It is noteworthy that every AIS for these 9 graphs with no symmetry contains an eigenvector of L 
with an integer eigenvalue, and every eigenvector of L with an integer eigenvalue is contained in 
an AIS. For example, (1, 1, 1, 1, 1, 1) G A c is an eigenvector of L with eigenvalue for any graph, 
and (2,2,— 1,-1,— 1,-1) G Ai is an eigenvector of L with eigenvalue 3 for the two graphs shown 
in Figure O 

Figure [5] also shows the bifurcation diagram for one of the graphs that has a non-constant AIS. 
The bifurcation diagram for the other graph is similar. We observe secondary bifurcations on the 
two primary branches bifurcating at the integer eigenvalues and 3. The secondary branch born 
of the constant solution at s = —3/2 contains solutions that are in A2, and these solutions have 
tertiary bifurcations to non- anomalous solutions. 

Finally, in Figure [5] one sees the phenomena of branch grouping by MI as s — > —00. We have 
observed this "grouping by MI" in all our experiments which use the schematic y(u) = \\u\\i, a 
superlinear /, and sufficiently negative s. The reason for the grouping is largely explained in |14j . 
where the asymptotic form of solutions in this realm takes on values in {0, c s , —c s } at each vertex, 
where f s (cs) = 0. The MI is computed by counting the number of nonzero components in the 
solution vector u, which also directly accounts for the value y{u) = 
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Figure 6. Cayley graphs of S3 (Example 17. 4|) . The graph on the top left is 
the Cayley color digraph Cayr/ 12 \ (23)}*§3- The graph on the right is a decorated 
Cayley graph with B3 symmetry. The highlighted vertices of the decorated graph 
correspond to the vertices of the Cayley color digraph. The bottom pictures show 
how the colored directed edges are replaced with decorated undirected edges. Since 
the generators are involutions, a pair of directed edges can be replaced by a single 
edge whose decoration does not encode the edge direction. 

7.4. A decorated Cayley graph of the symmetric group S3. Cayley graphs provide a way 
for us to generate a graph with a particular symmetry group. In Figure [6] we show how the Cayley 
color digraph Cay in 2) J (23)}§3 is used to generate a decorated Cayley graph with D3 symmetry. 
The symmetries and symmetry types for this graph are shown in Table [2j In Figure [7] we show un- 
condensed and condensed bifurcation digraphs containing all arrow types and nontrivial conjugacy 
classes. The uncondensed diagram has been annotated with contour plots to give visual cues as to 
the corresponding symmetries. The layouts and contour plots were all automatically generated by 
our suite of programs. 

The matrix L for this graph has the triple eigenvalue A4 = A5 = Xq = 3, whereas the irreducible 
representations of S3 are one or two-dimensional. Thus, the bifurcation point (u, s) = (0, 3) has an 
accidental degeneracy of Type 2, since the critical eigenspace E is the direct sum of two irreducible 
spaces (see Definition 15. ip . Such accidental degeneracy is a common feature of our experiments, 
since the matrix entries of the graph Laplacians are integers. The bifurcation of the constant branch 
at s = —1.5, seen in Figure HJ also has an accidental degeneracy of Type 2. 

7.5. The Cayley graph of Z5. We constructed a decorated Cayley graph of Z5 = (a | a 5 = 1) 
with 15 vertices. We conjecture that this is the smallest graph G with Aut(G) = Z5. This example 
is interesting for two reasons: First, there are non-EBL bifurcations with Z5 and Z5 x Z2 = Z10 
symmetry. Secondly, there are two inequivalent 2-dimensional irreducible representations of Z5 
with trivial kernels. 

Figure shows the irreducible representations of Z10 over M. All ten of the irreducible repre- 
sentations of Z10 over C are one dimensional: two are real, and eight are complex. Section 12.51 
describes how to construct the irreducible representations of Z10 over R. The three irreducible 
representations of Z5 over M. can be obtained by restricting to Z5. As described for general 
groups in Section [231 the irreducible representations of Ziq = Tq with a^ k \— 1) = —I are listed 

first. We choose this ordering so that V = ®^ = i and the isotypic components with k = 4, 5, 
or 6 satisfy vff = {0}. 
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Figure 7. Bifurcation digraphs for a decorated Cayley graph of S3 (Example IT.4[) . 
The digraph on the left is not condensed while the digraph on the right is condensed. 
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Table 2. Symmetries for a decorated Cayley graph of S3. The first column shows 
the isomorphism class of the elements in a condensation class; the second and third 
columns give the symmetry type and the symmetry. 



It is illuminating to describe the functions in Pp . In our layout of the graph, the generator a of 
Z5 acts as a rotation by 2tv/5. The vertices lie on three concentric circles, labeled by j G {1, 2, 3}. 

Let 9i be the angle of vertex Vi in the layout. Then = {u G V \ u% = Cj when Vi lies on circle j}. 

(2) (3) 
Similarly, Vp Q = {u G V \ Ui = Cj cos(6i) + dj sin(^j) when Vi lies on circle j} and Vp = {u G V \ 

m = Cj cos(26i) + dj sin(2#j) when Vi lies on circle j}. The dimensions of these isotypic components 

satisfy 3 + 6 + 6 = 15. 

There are exactly three symmetries in Q: Tq = (a, — 1) = Z5 x Z 2 = Z10, Ti = (a) = Z5, and 
r 2 = {1} = Zi. The symmetry types are singletons: S{ = {Ti} for i G {0,1,2}. The bifurcation 
digraph has exactly three arrows, a solid arrow So — > S± with label Z 2 , a dotted arrow 5o — > 5 2 
with label Z10, and a dotted arrow Si — >• 5 2 with label Z5. 

While there are three arrows in the bifurcation digraph, there are five bifurcation arrows for this 
graph, as shown in Figure El For example, the two bifurcation arrows To — )• T 2 , with labels k = 2 

and k = 3 both correspond to the single dotted arrow from So —> S" 2 with label Z10 in the bifurcation 

(2) 

digraph. The continuation solver needs to know if the critical eigenspace of the origin is in 

(3) 

or Vp , but from a theoretical point of view there is a bifurcation with Z10 symmetry in either 

case. The different nodal structures of the daughters bifurcating from the trivial solution shown in 

C^l (2) 

Figure |8| are explained by the above descriptions of Vp . When E C Vp , daughters change sign 
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Figure 8. The irreducible representations of Zio over M, the bifurcation digraph, 
and all five bifurcation arrows for the Cayley graph of Z5 (Example l7.5j) . The matrix 
of the rotation of M. 2 about the origin by 6 is denoted by R(9). The bifurcation arrows 
are labeled with the irreducible representation k. The arrow type and the group 
Ti/T' ik of the bifurcation are included in the arrow label to facilitate comparison 
with the bifurcation digraph. 



once, whereas when E C , daughters change sign twice. Similarly, the perturbations of the 
constant solution in the bifurcation with Z5 symmetry change sign once and twice, respectively. 

In Figure [9] we show bifurcations corresponding to the five bifurcation arrows. For a gradient 
system, a nondegenerate bifurcation with Zio symmetry creates 20 daughter branches in two group 
orbits of size 10, while a bifurcation with Z5 symmetry creates 10 daughter branches in two group 
orbits of size five. These bifurcations are similar to bifurcations with D10 and B5 symmetry, 
respectively. For example, a calculation shows that the normal form [12] for a gradient bifurcation 
with Z5 symmetry is g : C — > C, g(z) = Xz ± z\z\ 2 + (a + ib)z 4 , where A, a and b are real. The 
normal form for a bifurcation with D5 symmetry is the same, but with b = so that the real z 
axis contains solutions to g(z) = 0. The search directions in E that lead to daughter solutions are 
determined by the symmetry for D5, but remain a matter of trial and error for Z5. The situation 
is similar for the group Zio, except that z 4 is replaced by z 9 in the normal form. 

Our continuation solver had no trouble finding the bifurcating branches in this example. To test 
this, we ran the program that produced Figure [9] with several different seeds in the random number 
generator used in Algorithm [3] In most cases, all branches would have been found with / nc (2) = 5. 
In all cases, the default / nc (2) = 21 was large enough to find all of the branches. 

The keen observer will notice that the Z5 bifurcation with k = 2 in Figure [9] appears to have 
one branch bifurcating to the left and one branch bifurcating to the right, in contradiction to the 
fact that nondegenerate Z5 bifurcations have both of the daughter branches bifurcating the same 
direction. However, an extreme blow-up of the figure shows that both branches bifurcate to the 
left, but one branch has a fold point very close to the bifurcation point in addition to the visible 
fold point. 



7.6. The Cayley graph of the quaternion group Q. We are interested in the quaternion group 
Q = k I i 2 = j 2 = k 2 = ijk = —1) = {±1, ±i, ±j, ±/c} because it is the smallest group with 
a quaternionic representation (see Section 12.51 ) We find several examples of bifurcation with Q 
symmetry, which are interesting and complicated. Figure [10] shows the Cayley color graph of Q 
with generating set {i,j} and the corresponding decorated Cayley graph. 

Figure [TT1 shows the condensed bifurcation digraph computed by GAP. There are 14 symmetries 
in Q, and each of the symmetry types is a singleton: S% = {Ti} for all i. Since —1 G Q, we need to 
use the notation (g, ±1), with g G Q to denote elements of Q x Z2. The symmetries isomorphic to 
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Figure 9. Bifurcation diagram showing some of the solutions to Equation (UJ for 
the Cayley graph of Z5. The primary branches created at the first three distinct 
eigenvalues of L are shown, along with the daughter branches from two bifurca- 
tions with Z5 symmetry on the constant branch. The arrows and the black dots 
correspond to the six solutions shown in Figure [8l 




9 ■> gi 9 • — gi 



Figure 10. The Cayley graphs of Q (Example 17.6ft . The graph on the top left 
is the Cayley color digraph Cay^j^Q. The graph on the top right is a decorated 
Cayley graph, which has 48 vertices and 72 edges. The eight highlighted vertices 
in the decorated Cayley graph corresponding to the group elements g E Q. The 
bottom pictures show how the colored, directed edges are replaced with undirected 
edges. 
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Figure 11. The condensed bifurcation digraph for the Cayley graph of Q. Since 
the symmetry types are singletons, it is quite easy to deduce the bifurcation arrows 
from this figure using the description of the symmetries in terms of generators. For 
example, there are four bifurcation arrows emanating from r 2 : r 2 — > Tq, r 2 — > Tj, 
F2 r 10 , and T 2 -)• r !3 . 



Q are 

ri = ((i, 1), (j, 1)), r 2 = ((», -1), (j, 1)), r 3 = ((», 1), (j, -1)), and r 4 = ((*, -1), (j, -1)). 

These four subgroups of Tq = Q x Z 2 are not conjugate, but they are related by outer automorphisms 
and the four symmetry types are in the same condensation class, as seen in Figure [TTJ The six 
symmetries isomorphic to Z4 are 

r 5 = ((i, 1)), r 6 = ((*, -1)), r 7 = ((j, 1)), r 8 = ((j, -1)), r 9 = {{ij, 1)), and r 10 = {{ij, -1)). 

There are two symmetries isomorphic to Z 2 , namely 

r u = ((-1,-1)), andr 12 = ((-l,l)). 
The symmetry types Six and S12 are in different condensation classes. Finally, the trivial symmetry 

is ri 3 . 

Bifurcations with Q symmetry have a four-dimensional critical eigenspace E. These bifurcations 
are "highly non-EBL" , since all points in the critical eigenspace, except for the origin, have the same 
symmetry. Our continuation solver found examples of each of the 5 bifurcations with Q symmetry 
implied by Figure [TT] Tq — > Tn and Tj — > T13 for i = 1, 2, 3, and 4. This is the first time, to the 
best of our knowledge, that bifurcations with Q symmetry have been observed. Figure [12] shows 
the mother and one of the daughter solutions for a bifurcation T 2 — > T13. In a non-gradient system 
of differential equations, one would expect the bifurcation to create periodic solutions, but in our 
gradient system the daughters are solutions of V J s = that come in group orbits of size 8. 

We find the daughters by trial and error with repeated calls to the cGNGA function, as described 
in Algorithm [3j The number of random guesses can be modified by the user if index theory of 
something else suggests that not all daughter branches are found. In particular, the Poincare-Hopf 
index theorem [2] implies that 

(g) ^2 (-l) MI ("i s *- £ ) = ^ (_l) MI ( n ' s *+ £ ) ) 

(u,s*-e)eX (u,s*+e)eX 

where e is chosen so that the only bifurcation with s G [s* — £, s* + e] is at s = s*. In practice, we 
only need to sum over the mother and daughter solutions of the bifurcation. This invariant proved 
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MI 2 for s 6 (a*, A 2 ), MI 6 for s < s* MI 2 for s < s* 



Figure 12. Contour maps of solutions on the decorated Cayley graph of Q. The 
solutions are the mother (left) and one of the 80 daughters (right) of a bifurcation 
T2 — > I\3 with Q symmetry at s* ~ 0.328. The mother solution lies on the primary 
branch born at s = A2 ~ 0.347. The symmetry of the solutions can be determined 
from the value of u at the 8 highlighted vertices. For the mother solution, \u\ is the 
same at all highlighted vertices, and the action of i, shown in Figure [TO] switches 
the signs of u. Similarly, the signs are unchanged under the action of j. Hence the 
symmetry of the mother is T2 = ((i, — 1), (j, 1)). The daughter has trivial symmetry. 

particularly useful at bifurcations with Q symmetry. The MI on the mother branch changes by 
four, and each daughter has a group orbit of size 8, so in fact we sum Equation © over the set 
of nonconjugate daughters found by our continuation solver. The normal form for this bifurcation 
has not been computed, to the best of our knowledge. We have no theory predicting exactly how 
many daughters there are, or where in the 4-dimensional critical eigenspace the daughter solutions 
lie. ' 

Bifurcations with Q symmetry are extremely complicated. The daughters all go to the left in 
some examples, but some go to the left and some go to the right in others. We observed bifurcations 
with eight and ten nonconjugate daughter branches, for a total of 64 and 80 daughters, respectively. 
Some solutions had a very small basin of attraction in the cylinder, necessitating a large number 
of random guesses. The index invariant in Equation ([8]) was key in recognizing that some solutions 
were initially missing. However, index theory can never prove that all of the daughters have been 
found. 

One such bifurcation with Q symmetry occurred on the CCN branch. We found 10 nonconjugate 
daughters, all branching to the left, with MI 2, 3, 3, 3, 4, 4, 4, 4, 5, and 5, respectively. We had to 
increase num_no_changes = / n c(4) to at least 6,500 to satisfy index theory; an additional 360,000 
calls to cGNGA did not find any more daughters. For brevity, only the bifurcating CCN solution is 
shown in Figure [T2l 

7.7. The Petersen graph. We considered the well known Petersen graph in our experiments. 
The second eigenvalue for this graph Laplacian is of multiplicity 5 and there are 210 symmetries 
(20 types) of possible solutions to look for, which presents a challenge for our code. We are fairly 
confident that we have accurately followed a representative from each equivalence class of primary 
branches, and most if not all connected secondary branches. Under greater magnification, Figure [T3l 
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Figure 13. Full bifurcation diagram for the Petersen graph (Example I7.7p . There 
are high dimensional critical eigenspaces atu = 0,s = A2 = -- - = A6 = 2 and 
s = Xj = • • • = Aio = 5, as well at some secondary bifurcation points. Our code 
takes advantage of symmetry to search lower dimensional subspaces in order to 
efficiently find most solutions. The green MI 2 and blue MI 3 branches bifurcating 
from the second eigenvalue (see inset) contain a disconnected set of CCN solutions 
and are featured in Figure HH 

reveals that roughly 1300 points were used in following 75 branches, connected via 52 bifurcation 
points, with 3 reported (but not visually apparent) branch following failures. On our 3 GHz Linux 
workstation it took about 3 seconds to perform 1726 calls to tGNGA, with 4427 iterations at 2.565 
iterations per call; 433 cGNGA calls, with 2541 iterations at 5.868 iterations per call; and 83 secant 
calls, with 572 iterations at 6.892 iterations per call. 

A particularly interesting feature of the bifurcation diagram regards solutions with the minimum 
J value among all sign-changing solutions for that s parameter value, which are necessarily of MI 2 
(see [6]; for convenience we will call these CCN solutions here). In [7], we proved that CCN solutions 
exist up to A2, and in [17] we applied the GNGA to graphs and extended the CCN existence theorem 
(and related theorems) to graphs. We have since conjectured that there should exist a connected 
branch of CCN solutions for s € (—00, A2), but this numerical experiment indicates that this is not 
true. In Figure [U] we show the symmetry type £5 and Su primary branches bifurcating from the 
multiplicity 5 second eigenvalue. To the right of s* ~ 0.694, CCN solutions lie on the upper branch 
and have symmetry type S5, whereas to the left they lie on the lower branch and have symmetry 
type Su. In Figure [15] we provide contour plots of these two solutions at the crossover point s*, 
where both are global minimizers of J over the set of nontrivial sign-changing solutions. 

7.8. Dodecahedron. The space of functions on the dodecahedron graph admits 383 symmetries 
and 39 symmetry types. The default layout found by our program is close to the orthogonal 
projection of the 3-dimensional dodecahedron onto a plane parallel to a face (see Figure PT6|h 
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Figure 14. The CCN solutions for the Petersen graph, indicated by the thicker 
green lines, are not connected (see the inset in Figure [T3j) . To the right of s* ~ 0.694, 
CCN solutions have symmetry type S$ and lie on the upper branch, whereas to the 
left they lie on the lower branch and have symmetry type S±i. The CCN solutions are 
global minimizers of J over the set of sign-changing solutions, and always have MI 
2. This is a numerical counterexample to our previous conjecture that a continuous 
branch of CCN solutions exists for s < A2. 




Figure 15. Contour plots of two simultaneous CCN solutions and a MI 3 solution 
from the connecting branch for the Petersen graph, corresponding to the dots in 
Figure [TH The £5 solution has 12 elements in its symmetry group, six of which 
are visible in this layout. The Su and Sis solutions have symmetry groups of size 
four and two, respectively, although no nontrivial symmetries are visible in this 
layout. This layout, as well as the traditional layout of the Petersen graph were 
found automatically by our layout program. The layout used here makes more of 
the symmetries of the S5 solution visible. 

The dodecahedron features a Type 3 accidental degeneracy (see Definition I5.ip , which are rare 
in the examples we studied. At this bifurcation point the critical eigenspace is a direct sum of two 
1-dimensional irreducible subspaces lying in the same isotypic component. That is, K = {k} is a 

(k) \ (k) 

singleton set and dim(Vp\ D E) = 2 whereas dp = 1. 

The accidental degeneracy in this example can be explained using AIS. There is an AIS A a for 
the dodecahedron comprising of functions with Ui = Uj if vertex i and j are antipodal. It is well- 
known that the Petersen graph is the dodecahedron with antipodal points identified. Therefore, 
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the nonlinear operator VJ restricted to A a is the same as V J acting on functions on the Petersen 
graph. 

The operator V J\a u 1S equivariant under permutations on antipodal pairs of vertices. There are 
120 such permutations, since that is the size of the automorphism group of the Petersen graph. 
However, only 60 of these permutations are symmetries of the dodecahedron. We say that the other 
60 permutations are anomalous symmetries of VJ|^ a . 

The numerical results for the Petersen graph can be used to understand solutions for the dodec- 
ahedron in A a - At s* ~ 0.8727 there is a bifurcation with Z4 symmetry for the Petersen graph, at 
the same s value, there is a degenerate bifurcation of Type 3 in the dodecahedron for which the 
mother and daughters all lie in A a - One might expect that the symmetry of the mother in the 
dodecahedron would be isomorphic to Z4 x Z2, but in fact it is Z2 x Z2. Similarly, the symmetry of 
the bifurcation in the dodecahedron is not Z4, but rather there are two simultaneous bifurcations 
with Z2 symmetry. 

At the non-EBL bifurcation with Z4 symmetry in the Petersen graph, 8 daughters are created 
in two group orbits of size 4. One element in each group orbit is shown in Figure [TBI In the 
dodecahedron, the program also finds 8 daughters, but they lie in 4 group orbits of size two. We 
only show two daughters, since the others are conjugate under anomalous symmetries. 

7.9. Truncated icosahedron (soccer ball). We include a final example with more vertices and 
a large number of high multiplicity eigenvalues. The truncated icosahedron, made famous via the 
Buckminsterfullerene molecule, has 60 vertices. In Figure[T7]we display contour plots for 3 solutions 
from the second and third primary branches, near their bifurcations from the trivial solution at 
s = A2 = A3 = A4 ~ 0.2434 and s = A5 = • • • = A9 ~ 0.6972, respectively. This layout makes nearly 
all the symmetries of these solutions visible. It was found by our graph layout code, although it is 
not the layout with least complexity. 

8. Future Research Directions. 

Our suite of programs in their current state works well. We achieved our goal of taking an edgelist 
as input and automatically generating the symmetry information and solution data for PdE (jTJ) 
on the graph with that edgelist. The figures in this paper required only minor formatting of the 
raw results. We have tested our code on many other examples with a high degree of success. We 
successfully automated, for general graphs, the symmetry analysis found in our nonlinear snowflake 
code |l9j. The results encoded in the bifurcation digraph were used to follow most if not all 
bifurcating branches in these gradient systems. 

Bifurcation theory predicts that certain daughter branches must exist, but does not rule out the 
existence of other branches. Our continuation solver finds both types of daughter branches. In 
Example 17.61 we used index theory as an indicator that not all of the daughters had yet been found. 
It is an open problem to find more topological and variational theory to predict in general how 
many branches bifurcate, and where they lie in the critical eigenspace. Another open problem is 
to find a general theory of anomaly breaking bifurcations. 

Our focus in the current project was on large groups, not large graphs. For expedience we did 
not take advantage of all the methods used in |19J to speed up the calculations. That code was very 
efficient, using hardcoded symmetry shortcuts to reduce the number of integral calculations needed 
to define the linear systems required by Newton's method. The Hessian matrix h s (u) is a block 
matrix, with the number of zero blocks depending on the symmetry of u. In the snowflake code these 
zero blocks were not computed, but in the current work it was tolerated to perform the calculation 
of each element of h. We will implement this and other shortcuts for solving PdE on significantly 
larger graphs, specifically those obtained by discretizing PDE and using finite differences. Our next 
project will start with the application of our automated branch following algorithms to PDE on 
the square, which we first studied in [20J. 
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Figure 16. Contour plots of solutions to Equation (JTJ) for the Petersen graph (top 
row) and dodecahedron (bottom row) near the bifurcation point s* ~ 0.8727. The 
bifurcation has Z4 symmetry for the Petersen graph, but is a degenerate bifurca- 
tion of Type 3 for the dodecahedron. The MI information is valid on an interval 
(s* —e, s* + e) local to the bifurcation. The contour plots are obtained at s = 1.5 for 
the mother solutions (left column) and at s = — 2 for the daughter solutions. This 
layout of the Petersen graph shows the Z4 symmetry of the mother solution. The 
coordinates of the vertices in this layout were typed into a file rather than automat- 
ically generated by our layout program. Note that the four dots on the edge of the 
square are all the same size for the mother, but the Z4 symmetry is broken for the 
two daughter solutions. The bottom row shows corresponding solutions in the AIS 
of antipodal solutions A a on the dodecahedron. An interactive version of this figure, 
including three dimesional layouts of the Petersen graph and the dodecahedron, can 
be found at |http : //NAU . edu/ Jim . Swif t/PdE 



There are many PDE that merit an application of our PdE code. One area of interest is PDE 
on fractal regions. We propose to automatically generate large but finite pre- fractal graphs that 
in the limit converge to a fractal. Analytical and numerical research into the linear version of this 
problem has been done (see for example papers by R. S. Strichartz and A. Teplyaev and references 
therein), but nonlinear research where bifurcation is considered is absent from the literature. One 
could also investigate large graphs embedded in 2D manifolds such as the torus and sphere or 3D 
regions such as the cube. Generalizing L to approximate the Laplace-Beltrami operator is a related 
idea for future investigation. We are also interested in systems of PDE, and have made initial 
demonstration programs for computing solutions to several systems. Our code should perform well 
in investigating these types of problems, depending in part on our understanding of the underlying 
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Figure 17. Contour plots of three solutions to Equation ([T|) for the truncated 
icosahedron. The MI 2 CCN solution (left) has 20 symmetries, all visible in this 
layout. The front hemisphere is positive and the back is negative. The MI 5 solution 
(center) has 20 of 20 visible symmetries as well. It has a negative equatorial band 
separating front and back positive caps. The MI 6 solution (right) has 8 symmetries 
of which only 4 are visible. The nodal structure with two positive and two negative 
components is clear. All three solutions are very close to eigenvectors of L, which in 
turn resemble eigenfunctions of the PDE Laplacian —A on the sphere. Specifically, 
they have similar nodal structures as the spherical harmonics Y^o, an d Re (^2,2)1 
respectively [5]. 

theory for systems. Another area of future research is to borrow from the established linear graph 
theory and our ghostpoint ideas from |18j to accurately enforce alternate boundary conditions to 
our PdE (and hence PDE) code. 

We are interested in the existence, multiplicity, and nodal structure of solutions to nonlinear 
elliptic PDE. Thus, we will continue to perform experiments to support conjectures in the analyt- 
ical theory for PDE, seeking a better understanding of the underlying variational structure. We 
found several interesting phenomena in the PdE examples and seek to determine if they persist for 
PDE. For example, in Section 17.71 one sees an example where the CCN branches of solutions for a 
PdE are disconnected, contrary to our conjecture that a continuum of such solutions exists for all 
s < A2 (assuming the standard subcritical/superlinear hypothesis found in [1, 6j). Furthermore, 
in Example 17.21 we found strong numerical evidence that not all symmetry types are present in 
the solution set X for PdE. It would be instructive to find examples of PDE with similar features. 
Finally, we would like to find a PDE result analogous to the grouping by MI property commented 
on in Example 17.31 
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